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Abstract. We study a two-dimensional kinetic Ising model with Swendsen-Wang 
dynamics, replacing the usual percolation on top of Ising clusters by explosive 
percolation. The model exhibits a reversible first-order phase transition with hysteresis. 
Surprisingly, at one of the transition flanks the global bond density seems to be equal 
to the percolation threshold. 
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Three years ago Achlioptas et al. discovered that a slight modification of standard 
percolation on a complete graph turns the continuous transition into an abrupt one [I] . 
This surprising phenomenon, marketed as "explosive percolation", attracted a lot 
of attention and was studied on scale-free networks (2j|3), regular two-dimensional 
lattices [4j, and in various other situations [5 12 . Initially the explosive transition was 



believed to be discontinuous, but later studies revealed that the transition is extremely 



sudden but actually continuous 13 -17 



While in ordinary bond percolation each bond is set independently with 
probability p, the Achlioptas process uses a dynamical procedure, where the bonds 
are placed in sequence. To this end one randomly selects two pairs of neighboring 
lattice sites which are not yet connected by a bond. Depending on the actual bond 
configuration, a weight is assigned to each of the pairs by multiplying the sizes of the 
clusters to which the two lattice sites belong (or squaring the size if both sites belong to 
the same cluster). If the weights of the two pairs are different, the pair with the lower 
weight is connected by a new bond. Otherwise, if the weights are equal, one of the pairs 
is randomly chosen and connected. This process is repeated until the density of bonds 
exceeds the value of the control parameter p. 



Explosive Ising 



2 



rPr tJ n 

i — i i • i — { {—i 

n 



n 



^4 

♦ ft 



Figure 1. Left: Dynamics of explosive bond percolation on a square lattice. 
Two vacant bonds A, B (shown as dashed lines) are randomly selected. As bond 
A would connect two clusters of sizes 2 and 4, its weight is 2 • 4 = 8, whereas 
bond B has the weight 3-3 = 9. The bond with the lower weight is set and 
the procedure is repeated until the density of bonds exceeds the parameter p. 
Right: "Explosive" Ising model with two spin domains (red and green). Here 
a Swendsen-Wang update uses the same dynamical rule for placing bonds but 
only within domains of equally oriented spins. 



As the dynamic rule is denned in such a way that new connections are preferentially 
added between small clusters, it is plausible that the transition point is shifted to 
higher values of the percolation probability. For example, on a two-dimensional 
square lattice, where the transition of ordinary bond percolation takes place at 
pperc _ the Achlioptas-dynamics described above shifts the transition point to 
^expi _ 0.526562(3) [14]. Surprisingly, the transition is not only shifted, but it changes 
also qualitatively, exhibiting a sudden transition, where a large percolation cluster is 
formed. This is why the transition is called "explosive" . 

In this Letter we demonstrate that the concept of explosive cluster growth can 
also be applied successfully to kinetic spin models with cluster dynamics. Here we 
will focus on a particularly simple case, namely, the two-dimensional Ising model with 
Swendsen-Wang dynamics |18|. As will be shown below, the use of an explosive cluster 
dynamics turns the continuous Ising transition in a discontinuous one, preserving the 
Z 2 -symmetry of the model. As this modification of the dynamics is expected to break 
detailed balance, the stationary state of such an explosive Ising model will no longer be 
an equilibrium state. 

Definition of the model: 

Before defining the explosive Ising model let us briefly recall the Swendsen-Wang (SW) 
algorithm applied to the Ising model [18]. The SW cluster algorithm is a dynamical 
update rule which works as follows. At first percolation clusters are grown within the 
actual spin clusters. For each of these clusters one chooses a random number ±1 which 
is then assigned to all its spins. More specifically, in case of the Ising model the SW 
update for a given configuration of spins consists of the following steps: 
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(i) Remove all bonds. 

(ii) Set all bonds between equally oriented nearest-neighbor spins with probability p. 

(iii) Identify all clusters of sites which are connected by bonds. 

(iv) For each of these clusters generate a random number ±1 with equal probability and 
assign it to all its spins. 

It was proven that this dynamics evolves into the equilibrium state of the Ising model 
without critical-slowing-down and that the percolation probability is related to the 
temperature by p = 1 — e -2J / fcsT . 

Let us now modify the SW dynamics of the Ising model by replacing ordinary 
with explosive percolation. Since bonds can no longer be distributed independently, one 
first has to determine the number of potential bonds n. This results into the following 
"explosive" Swendsen-Wang dynamics for the Ising model (see right panel of Fig. [T]): 

(i) Clear all bonds. 

(ii) Count the number n of links between neighboring spins of the same orientation. 

(iii) Select two vacant bonds between equally oriented spins and determine their weight 
according to the product rule in the same way as in explosive percolation. 

(iv) Place a bond at the link with the lower weight (if the weights are equal select one 
of them randomly). 

(v) Repeat (iii) and (iv) until the number of bonds exceeds np, where p G [0, 1] is the 
control parameter of the model. 

(vi) Finally, as in ordinary SW dynamics, assign a random spin orientation to each of 
the clusters. 

Fig. [2] shows typical spin configurations of the ordinary Ising model compared with those 
of the explosive variant defined above, increasing the parameter p from 0.55 to 0.59. For 




Figure 2. Typical snapshots of spin configurations for the ordinary and the 
explosive Ising model with p increasing from 0.55 to 0.59. 
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Figure 3. Order parameter as a function of the percolation probability p for 
various models discussed in the text. From left to right: (a) Average density 
of the largest cluster in ordinary bond percolation on a 3000 x 3000 lattice and 
periodic boundary conditions. The dotted line indicates the transition point in 
the limit L — > oo. (b) Largest cluster size in explosive percolation as defined 
above on a periodic 500 x 500 lattice, (c) Hysteresis of the order parameter 
\M\ in the explosive Ising model on a square lattice with 512 x 512 sites and 
periodic boundary conditions, (d) Same data for the ordinary two-dimensional 
Ising model (analytical result). 



small values of p both models are in the disordered phase, although the domains in 
the explosive variant seem to be somewhat smaller. As p is increased, the Ising model 
displays a continuous transition at p^ mg = 2 — \/2 w 0.585786 into a partially ordered 
state, while the explosive variant switches suddenly to an almost completely ordered 
state at some value between 0.57 and 0.58 (see below). 

Order parameter hysteresis: 

In explosive percolation the process of adding new bonds is not reversible, i.e. in a single 
realization the density of bonds can only be increased. In the present model, however, 
the explosive cluster dynamics is part of a continually repeating update procedure, 
allowing one to increase and decrease the percolation probability while the model evolves 
in time. 

Since a spin model with cluster dynamics switches frequently from positive to 
negative magnetization in the ordered phase, the appropriate order parameter is the 
absolute value of the magnetization M. The solid lines in Fig. [3] show how \M\ varies 
with the control parameter p on a lattice with N = 512 x 512 sites and periodic boundary 
conditions. As can be seen, the transition seems to be discontinuous, exhibiting a 
pronounced hysteresis with the flanks located at the values p Cj i and p Cj2 (see Table 1). 



Explosive Ising 



10000 





30000 



10 



1 1 1 1 1 


I 


1 1 1 E 


; N=24x24 


— p= 


=0.555 a 




- P= 


=0.560 L 




- p= 


=0.565 / ! 




— P= 


=0.570 Ift 






J I 









0,2 



0,4 0,6 

IMI 



" 




2000 



N=L 



Figure 4. Left: Spontaneous flipping of the order parameter \M\ in a system 
with 24 x 24 sites inside the hysteresis loop at p = 0.563. Middle: Distribution 
of the order parameter averaged over long time for various values of p in a 
24 x 24 system. Right: The average time T(N) between two flips at p = 0.563 
grows exponentially with the system size. 



The hysteresis may be explained as follows. For p 0.5 the Achlioptas selection rule 
generates smaller clusters than in the usual Ising model. Since the percolation clusters 
are confined to domains of the same spin orientation, their growth is suppressed like 
in a finite system with open boundary conditions, explaining why nothing happens at 
the usual threshold p^ xpl ~ 0.527 of explosive percolation. However, as the size of the 
Ising domains grows with p, one eventually reaches a point p c ^ where this stabilization 
mechanism breaks down. Here one observes the immediate formation of a system- 
spanning cluster, self-sustained from inside by the as from now super-critical Achlioptas 
process. Lowering p again the system first remains in the ordered state until the explosive 
percolation transition is reached from above. As the percolation process takes place on 
a slightly porous support, the transition point is already reached at a value p c l which 
is somewhat higher than the usual transition point of explosive percolation p^ xpl . 

Spontaneous flipping and order parameter distribution: 

As usual in systems with a first-order phase transition, we find that the order parameter 
in finite systems flips occasionally between the disordered and the ordered phase (left 
panel of Fig. As expected, the corresponding probability distributions of the order 
parameter (see central panel) show two maxima separated by a valley, whose relative 
sizes depend on p. We also verified that the flipping time T(N) diverges exponentially 
with the number of sites N = L 2 , meaning that the ordered and the disordered phase 
are both thermodynamically stable in the coexistence region p Cj i < p < p c ^- 
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explosive 
percolation 
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P r c = 1/2 


pf? 1 = 0.526562(3) 


p Ci i = 0.5516(3) 


p c2 = 0.5725(3) 


pf n § = 0.585786 



Table 1. Transition points for various models discussed in the text. 
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Density of equally oriented neighbors at the transition flanks: 

Let b be the total number of bonds generated in the update procedure. We finally 
present the surprising conjecture that the global bond density 

' = 2* m 

is exactly equal to p^ xpl at the left transition flank in the thermodynamic limit. 

Remarkably, a similar relation holds at the transition point of the ordinary Ising model 
on a square lattice [19]. To see this let us consider the number n of equally oriented 
neighboring spins, which is related to the internal energy E = —JJ2(ij) s i s j °y 

E + E 

n = , 2 

2J ' K J 

where Eq = —2NJ the lowest possible energy where all spins are parallel. As bonds are 
randomly set between equally oriented spins with probability p, we have b ~ np so that 
the global bond density is given by 

q = M im OAT' ( 3 ) 

Since the internal energy of the Ising model on a square lattice at the critical point 
p isin g = 1 _ e -2j/k B T c = 2 - v/2 is known to be E c = -y^NJ, we find that 



2N 4 v ; 2 v ; 

This means that the transition in the ordinary Ising model with Swendsen-Wang 
dynamics takes place when the global bond density reaches the value 1/2 - the same 
value as the critical threshold for percolation without Ising spins. The actual transition 
point p* smg pa 0.5858 is of course larger since bonds between different spin clusters cannot 
be set so that an enhanced probability is needed to reach a global bond density of 1/2. 

Surprisingly, we find numerically that the same relation holds at the left transition 
flank of the explosive Ising model, where the bond density g Cj i reaches the value p^ xpl of 
explosive percolation without Ising spins (see Fig. [5]). However, at the right flank the 
transition occurs at a lower value q Ct 2 ~ 0.49 which may tend to 1/2 in the limit of very 
large lattices. 




Figure 5. Left: Global bond density q set during the SW update as a function of 
the percolation probability p in the ordinary Ising model (left) and the explosive 
variant (right). The horizontal dotted lines indicate where q = pc CTC = 1/2 and 
q = pT pl = 0.526562(3). 



Concluding remarks: 

In this Letter we introduced and studied an explosive variant of the Ising model, arriving 
at two conjectures: 

• Even though the Achlioptas-process was found to be continuous, we think that the 
transition in the explosive Ising model is genuinely discontinuous because it exhibits 
a very clear hysteresis. 

• The global density of bonds in the SW update at the left transition flank of the 
hysteresis is conjectured to be equal to the threshold p^ xpl of explosive percolation 
without Ising spins. An analogous relation holds already for the usual Ising model. 
At the right transition flank the bond density is lower, probably close to 1/2. 

It would be interesting to explore these conjectures more deeply and to apply similar 
ideas to other spin systems with cluster dynamics. 
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